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By using molecular dynamics simulation, formation mechanisms of amorphous carbon in particular sp 3 
rich structure was researched. The problem that reactive empirical bond order potential cannot represent 
amorphous carbon properly was cleared in the transition process from graphite to diamond by high pressure 
and the deposition process of amorphous carbon thin films. Moreover, the new potential model which is 
based on electron distribution simplified as a point charge was developed by using downfolding method. As 
a result, the molecular dynamics simulation with the new potential could demonstrate the transition from 
graphite to diamond at the pressure of 15 GPa corresponding to experiment and the deposition of sp 3 rich 
amorphous carbon. 



1. Introduction 

There is a lot of amorphous material in the world. In the past, almost all of them had been 
classified into 'amorphous' just because they do not have crystalline structures. However, ma- 
terials can be treated in nano-meter scale currently, and then the fact that the amorphous 
materials with different structures and properties elements in spite of same constituent atomic 
elements exist is emerged. The reason for the structural diversity of amorphous carbon mate- 
rials is that an sp 3 structure which is a metastable state is as much stable as an sp 2 structure 
in carbon. From this reason, one of important parameter to classify amorphous carbon ma- 
terials is the sp 3 /sp 2 ratio, which is the ratio of the number of sp 3 carbon atoms to that of 
sp 2 carbon atoms. Actually, in also the most well-known classification of amorphous carbon 
materials proposed by W. Jacob, W. Moller and Robertson, 1 ' 2 ) the sp 3 /sp 2 ratio was selected 
as a classification parameter with hydrogen content rate. However, the amorphous carbon 
materials which are same in sp 3 /sp 2 ratio and hydrogen and but differ in physical property 
such as strength exist. The classification by these two parameters is, therefore, not enough 
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to understand amorphous carbon materials. In other words, if more detail classification is 
obtained, it can lead to the discovery of the new amorphous carbon materials which have 
physical property functionality such as friction property, electric conduction property, semi- 
conductor property, and so on. Aiming to establish new classification of amorphous carbon 
materials, we have investigated by using molecular dynamics (MD) simulation. MD simula- 
tion is powerful method because it can demonstrate amorphous structures in atomic scale 
directly. Actually, our previous research cleared that the amorphous carbon materials created 
by deposition process or annealing process differ in the orientation of covalent bonds although 
they are similar in density, sp 3 /sp 2 ratio and radial distribution function. 3 ) It follows from this 
that to clear the difference of amorphous carbon materials, their formation processes should 
be understood simultaneously. 

The purpose of this study is to investigate the amorphous carbon materials of sp 3 rich 
structures. In general, almost of all carbon materials have sp 2 rich structures and diamond, 
which is typical material of an sp 3 rich structure, is generated in high pressure environment. 
However, diamond like carbons (DLCs) are created by deposition process in the condition of 
low pressure and then they have sp 3 rich structures. It is not clear how DLCs can have a lot 
of sp 3 carbon atoms in amorphous structure. 

In the previous work, we investigated amorphous carbon deposition by using MD simu- 
lation with the second generation reactive empirical bond order (REBO2002) potential. 4,5 ) 
As a result, relationship between deposition ratio and H/C ratio of incident atoms onto a 
surface, and then the relationship agrees with experimental result. 6 ) After that, using similar 
MD simulation we have tried to clear the difference between formation process of sp 2 and 
sp 3 rich amorphous carbon deposits. Simultaneously, we have tackled graphite to diamond 
transition under high pressure because the transition is regarded as a structure change from 
a pure sp 2 structure of graphite to a pure sp 3 structure of diamond via sp 2 rich and sp 3 rich 
amorphous carbons structures. It is important to study from various viewpoints to elucidate 
the complicate mechanism of amorphous carbon. However, we encountered the obstacles that 
sp 3 rich amorphous carbon cannot created in MD simulation and graphite cannot change into 
diamond by more high pressure than 100 GPa, where experimental transition pressure was 
reported by about 15 GPa. 7 ) We consider the reason of this problems is that the empirical 
function forms in the REBO2002 potential have trouble representing potential energy for 
amorphous structures. 

Here, we note the empirical function forms in the REBO2002 potential. In general, almost 
all model of potential in MD of real atomic system, in which dominant interaction is covalent 



bonding especially, is classified into 'bond order (BO)' potential. BO potential is composed 
of the functions of the length of covalent bond, the bond angle which is the angle between 
adjacent two covalent bonds, the dihedral angle which is the angle between adjacent two planes 
formed by continuous three covalent bonds, and so on. As a extension for BO potential to treat 
new connecting and cutting of covalent bonds, Brenner et al. added the term of functions of 
the number of surrounding atoms within a cutoff length. This surrounding term represents the 
energy difference among molecular structures empirically, and then the potential energy of a 
reaction intermediate molecule was given by interpolation of energies of molecules before and 
after reaction, which is called reactive empirical bond order (REBO) potential. 8 ) The REBO 
potential is most often employed in the study of carbon nano-materials since the discovery 
of carbon nanotube and fullerene. Moreover, the second generation model of REBO potential 
(REBO2002), which had been improved over a period of twelve years by Brenner et al., 4 ) can 
calculate not only reasonable energies for almost all carbon and hydro-carbon molecules. 

Thus, the REBO2002 potential well represents stable molecule and lattice structures in 
potential energy. However, meta-stable structures such as amorphous structures are hardly 
considered in the REBO2002 potential and general BO potentials. Therefore, in the present 
work, we propose new potential model for amorphous carbon material. The new potential 
model is based on consideration of electron distribution like molecular orbital, and is not BO 
and REBO potentials. Consequently, in the MD simulation with the new potential model, the 
creation of sp 3 rich amorphous carbon in a deposition process and the graphite to diamond 
transition under pressure at 15 GPa could succeed. 

2. Electron Order Potential 

2.1 Potential function form 

As a function form of the new potential, BO and REBO type functions are not employed. 
Then, from natural viewpoint, total potential energy [/tot is composed of nucleus- nucleus 
repulsive energy U nn , nucleus-electron attractive energy U nc , and electron-electron repulsive 
energy U cc : 



The nucleus-nucleus repulsive energy U nn has a function form of screened Coulomb inter- 
action given by 



U nn + U nc + U ( 
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where Yli j^i means summation in terms of pair of i-th and j-th nuclei, Qi and Qj are core 
charge of nuclei which are screened by core electrons, r^- = \ri — rj\ is a distance between the 
i-th and j-th. nuclei, and / nn (ni) is cutoff function. The detail of cutoff functions is described 
later. 

To construct the nucleus-electron attractive energy U ne and electron-electron repulsive 
energy U ee , electron distribution like molecular orbitals is considered. Of course, the electron 
distribution is simplified as a point charge located on the position of atom (nuclei) r» and the 
center position between pair atoms = (rj + rj)/2. When the i-th and j-th atoms are close, 
electric charge transfer from the position of the i-th atom to the center position between i-th 
and j-th atoms is defined by 

««(r«) = kC-^^) (3) 
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if r > 1, 

/q( r ) = \ (~6r 2 + 15r - 10)r 3 + 1 if 1 > r > 0, (4) 
else, 
where single electric charge is normalized as 1. In this definition, if > r cs + r e( j, qij{fij) = 0. 
That is, electron is unpaired electron completely. Therefore the electric charge of unpaired 
electron located on the position of the i-th atom is 

where qj is number of valence electrons depending on the atomic element of the i-th atom. 
The nucleus-electron attractive energy U ne is divided into four terms as follows: 

U nc = U s ne + U% + CC + U£ (6) 

The term U* e is the attractive energy between nucleus and the nearest shared electron pair 
defined by 

U ne = -Yl K ™ exp (-< e fa - n\) Qi [g y (r^) + (fcifai)] ■ (?) 

The term U^, is the attractive energy between nucleus and the second nearest shared electron 
pair defined by 

^ne = - K ™ eX P (~ a ne \^ik ~ Tj |) Qj fefc(^fc) + ffifeifai)] fnl( r ij), (8) 

where /nc( r «i) is cutoff function. The term U£ e is the attractive energy between nucleus and 



the nearest unpaired electron denned by 



(9) 



where the electric charge of unpaired electron is given by eq. (5). The nearest unpaired electron 
is located on same position to nucleus. The term U£ e is the attractive energy between nucleus 
and the second nearest unpaired electron defined by 
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where /^(fy) is cutoff function. Because the second nearest unpaired electron is located on 
same position to the second nearest nucleus, the distance between nucleus and the second 
nearest unpaired electron is same to that between pair of nuclei. 

The electron-electron repulsive energy U ee is modeled by using the electric charge of shared 
electron pair of eq. (3) and that of unpaired electron of eq. (5) similarly to U ne . The electron- 
electron repulsive energy U ee is composed of the following three terms: 
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The first term U£ e is repulsive energy between two electrons of shared pair defined by 

Ul = \j2Klq ij (r ij )q ji (r ji ) ) (12) 

The second term C/|| is repulsive energy between adjacent two shared electron pairs defined 
by 

U ee= Yl K el^ x P{-af e \fij - fik\) [Qij(rij) + tyifai)] [ftfcfofc) + q k i(r ki )] . (13) 
i,jjti,k>j 

The third term £7™ is repulsive energy between shared electron pair and the nearest unpaired 
electron defined by 



= A'cc exp (-a™ \fij - n\) fajirij) + 
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There are three cutoff function f nn , f%£, /~ in the above definitions. In the present work, 
all cutoff functions are unified by using eq. (4) as 



fnn(r) = f£(r) = / n u2 (r) = / q ( 
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Except for the nucleus- nucleus repulsive energy U nn in eq. (2), all terms have no Coulomb 
potential form. This reason is that U nc and U ee are regarded as interaction energy of not 
so much classical electron particle but electric charge distribution due to quantum electron 



orbitals in the present model. Although the potential model in the present work is composed 
of only the above terms, higher quality potential model can be created by adopting nucleus- 
electron attractive terms and electron-electron repulsive terms between more distant positions. 

2.2 Downfolding method 

Ideal potential model should be able to calculate proper inter-atomic interaction energy in 
arbitrary atomic geometry. In other word, the potential model which can calculate inter- 
atomic interaction energy in many kinds of atomic geometries with accuracy is regarded as 
good model. Currently, development of the quantum chemistry which is represented by density 
functional theory (DFT) are made it possible to calculate the inter-atomic interaction energy 
with a high degree of accuracy. Then, by converging the parameters in potential function so 
as to reduce the difference between inter-atomic interaction energies calculated by potential 
model and DFT in many atomic geometries, we can obtain good potential model. Yoshimoto 
proposed the downfolding method to develop potential model from these point of view. 9 ) In 
the present work, the parameters in the function form given in §2.1 are optimized by the 
downfolding method as follows. 

Potential is, of course, a function of atomic geometry {r} = {r\,V2, ■ ■ ■}. If parameters in 
potential function are regarded as variables, the potential is also a function of the parameters. 
The parameters in potential function is described as {a} = {cq, a2> ■ ■ ■}) and the potential 
function is described as U({r},{a}). The variance of difference between inter-atomic inter- 
action energy calculated by using potential model and DFT in terms of iV kinds of sample 
atomic geometries is given by 

* = jf E P({ r h > W) " E rei {{rh)} 2 (16) 
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where {r} i means the i-th sample atomic geometry and the reference energy £? re f({r}J is the 
inter- atomic interaction energy calculated by using DFT. As so to reduce 3>, the parameters 
{a} is optimized. 

One way to optimize the parameters {a} is to solve the following evolution equation on a 
virtual time t: 

§-^- 2 ^ ^^W^ t c,<Mi ■{•*»- *-«->•>] ■ (i7 » 

where the parameters {a*} = {a\, o^, cdots} indicate the parameters at the virtual time t and 
the coefficient c is weight for time evolution. The parameters {a*} when t — > oo is optimal 
parameters. In the numerical solvent, this evolution equation is replaced with the difference 



equation given by 
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In actual operation, how to prepare the sample atomic geometries is important. Yoshimoto 
proposed that the sample atomic geometries are selected from animation snapshot in prelim- 
inarily performed MD simulation with the potential function C/({r},|o }) which employs 
temporary parameters {a }- For example, if the preliminary MD simulation is executed in 
canonical ensemble scheme, the sample atomic geometries are generated according to canon- 
ical distribution. Yoshimoto adopted multi-canonical ensemble scheme as a preliminary MD 
simulation. Moreover, the iteration of the above sequence of the downfolding method, replac- 
ing temporary parameters {a } with the optimal parameters given by the previous sequence, 
is an effective technique to obtain more optimal parameters. This iteration is also effective to 
prevent the parameters from being trapped in local minimum. 

In the present work, the optimized parameters {a} in the above sequence are correspond- 
ing to fifteen parameters i^nn? ^nnj -^-"ne? ^ne? ~^ne> ^ne> -^ne? ~^ne > ^ne> -^ee? -^-ee' ^ee> ^ee ' 

a™ } in the function form defined in §2.1. Sample atomic geometries were created by the three 
kinds of preliminary MD simulation with REBO2002 potential. One third of sample atomic 
geometries were selected from graphite structures vibrating at high temperature, the next one 
third of sample atomic geometries were selected from diamond structure similarly, and the 
other one third of sample atomic geometries were selected from amorphous structures made 
by pressing graphite in the preliminary MD simulation. The reference energy E^f^r}^ for 
the sample atomic geometries were calculated with 'Open source package for Material eX- 
plorer' (OpenMX), 10 ) which is numerical software based on DFT. The generalized gradient 
approximation (GGA) with Perdew-Burke-Ernzerhof (PBE) functional 11 ) was employed as 
exchange-correlation potential. Norm conserving pseudo-potentials 12 ) were used for approx- 
imation for core electrons and nuclei. Electron orbital was represented by numerical orbital 
with pseudo-atomic localized basis functions. 

As a result, the optimal parameters were obtained, where the iteration of the downfolding 
method was not performed. The optimal parameters and fixed parameters are shown in Table 
I. 



3. Application 

3.1 Transition from graphite to diamond 

The first results by the new potential was that the transition from graphite to diamond 
could be well imitated by MD simulation. In this MD simulation, it was investigated that the 
lattice structure changes from a graphite structure into a diamond structure as the pressure 
of system increases gradually. The system had a 24 carbon atoms, and initially they construct 
the graphite structure which was composed of three graphene layers stacked as a 'ABC form. 
The simulation box was under the periodic boundary condition. The pressure of system was 
controlled by the Andersen's method [ref], where a setting pressure increased from 3 GPa to 
15 GPa in the first 4 x 10 _12 s, and then it was kept at 15 GPa during the next 4 x 10~ 12 s. 
Though the size of simulation box was initially 0.492 x 0.426 x 1.004 nm 3 , it was changed 
according to the pressure of system on the basis of the procedure of the Andersen's method. 
Different point from the original Andersen's method was that the pressure and the size of 
simulation box were controlled independently for each x, y, z directions. The temperature 
of the system was also controlled by the Langevin thermostat method, in which a setting 
temperature was 1000 K and a friction coefficient is 1.0 x 10 1 4 s -1 . The time step of the MD 
simulation was 4.0 x 10~ 17 s. After 2 x 10 5 steps, the MD simulation finished. 

Figure 1 shows the change of pressures in the MD simulation. From this figure, it was 
confirmed that the pressure acting actually on carbon atoms well follows the setting pressure 
of the Andersen' s method. In this simulation, we could demonstrate the transition from 
graphite to diamond at the pressure of 15 GPa, which is agreement with experimental report, 7 ) 
as shown in Fig. 2. Quantitative representation of the transition from graphite to diamond 
was given by the change of ratio of the numbers of sp 2 and sp 3 carbon atoms to that of all 
carbon atoms as Fig. 3, which are here called sp 2 and sp 3 ratios, simply. As the pressure 
increased, the sp 2 ratio decreased and the sp 3 ratio increased. When the pressure reached 15 
GPa, the sp 3 ratio also achieve to become 1.0. That is, carbon atoms became a perfect sp 3 
material. From the simulation snapshot shown in Fig. 2(last), it was confirmed that when 
the sp 3 ratio was 1.0, the carbon atoms certainly constructed a diamond structure, not an 
amorphous structure. 

3.2 Deposition of amorphous carbon 

Using the present new potential model, we tried to simulate the deposition generating an sp 3 
rich amorphous carbon. The manner of the MD simulation of deposition was similar to our 
previous work. 6 ) The substrate was prepared by piling the amorphous carbon blocks which was 



created by the MD simulation for deposition on the REBO2002 potential. 3 ) The sp 3 ratio of 
the substrate was about 20 percent. The size of the surface of the substrate was 2.019 x 2.186 
nm 2 and then the substrate was put on the simulation box which follows periodic boundary 
condition in the parallel direction to the surface of the substrate. The carbon atoms initially 
putted on simulation box as the substrate were connected to the Langevin thermostat with a 
friction coefficient of 1.0 x 10 1 4 s _1 to control temperature, where the motion of carbon atoms 
located in the range of up to 0.12 nm from the bottom position was fixed during simulation. In 
one simulation, 1000 carbon atoms as a source were continuously injected from 3.4 nm above 
the initial surface of the substrate every 1.0 x 10 12 s. Injection position in parallel direction 
to the surface of the substrate was determined randomly under uniform distribution. These 
MD simulations were executed varying the injection energy from 1 eV to 200 eV. 

As a result of the MD simulation using the present new potential, we could obtain the sp 3 
rich amorphous carbons as deposits. Figure 4 shows sp 2 and sp 3 ratios of only the injected 
carbon atoms staying in the deposit, except for that constructing initial substrate and that 
reflected from a surface. From this figure, it was confirmed that the sp 3 ratio in the case of 
the present new potential was grater than that of the REBO2002 potential for whole injection 
energies. 

The penetration depth of the injected carbon atoms was also investigated as shown in Fig. 

5. Here, the penetration depth was defined as the mean of the distance between the surface 
position at the moment of each injection and the final position of each injected carbon atoms. 
In the case using the present new potential, the penetration depth increased as the injection 
energy increased higher than 20 eV, while it was independent of the injection energy for lower 
range. On the other hand, in the case using the REBO2002 potential, the injected carbon 
atoms hardly penetrated from the surface. The penetration depth in case of the injection 
energy of 200 eV is even corresponding to three times of thickness of diamond and graphite 
monolayer. 

Deposition rate was also estimated by the time evolution of the thickness of deposit as Fig. 

6. Here, the thickness was defined by the difference of the surface position from the minimum 
position while simulation. The reason why the thickness decreased for a while from start of 
simulation is by the contraction of substrate. The substrate composed of stacked amorphous 
block, which was created by using the REBO2002 potential, had an sp 2 rich structure and is 
unstable state for the present new potential, the thickness decreases. The simulation in the 
case of high energy injection used many amorphous blocks to prepare the substrate, and the 
decreases of thickness was larger. After the contraction of substrate, the thickness increased 



linearly. By this linear growth, the deposition rate could be evaluated as the gradient of the 
increase of the thickness. It is understand that the deposition rate was independent of the 
injection energy of less than 100 eV because sputtering hardly occurred. The case of the 
incident energy of 100 eV or more, we should calculate more long time. 

4. Discussion 

In the present paper, we proposed the new potential model. In fact, before developing the new 
potential model, we had tackled the modification of the REBO2002 potential to represent the 
sp 3 rich amorphous carbon. Especially, in terms of the three body terms depending on bond 
angles, many function forms were tried and parameters were optimized by using the down- 
folding method every time. However, we could not obtain well modified potential to achieve 
the transition from graphite to diamond due to high pressure. Compromising optimization 
in parameter, the graphite kept its structure even of pressure over 100 GPa or it changed to 
highly-oriented structures. From this experience, we consider that approximation of interac- 
tion of covalent bonding system with the bond angle, dihedral angle and so on, which are 
historically employed by bond order type potential models, does not give us well representa- 
tion for amorphous structures, which are metastable structures. For the REBO potential, it 
was reported that the elastic constant of diamond is not also properly evaluated. 13 ) 

On the other hand, the bond angle and the dihedral angle were not employed as the 
variables of the function form of the present new potential. Although simplified to point charge, 
electron distribution was considered and then potential function consists of nuclear-nuclear 
repulsion, nuclear-electron attraction and electron-electron repulsion according to that natural 
consideration. We carried out the downfolding for this function form, and we could smoothly 
optimize all parameter. Consequently, we could amazingly-easily achieved the transition from 
graphite to diamond by high pressure. In Fig. 7, potential energies calculated by the present 
new potential, the REBO2002 potential and DFT were compared in terms of several carbon 
structures in the MD simulation with the present new potential for the transition from graphite 
to diamond in §3.1. The potential energies by the new potential model are close to those 
by DFT, while those by REBO2002 were significantly different from those by DFT. This 
comparison implies the advantage of the present new model in treatment of amorphous carbon 
structures. 

Furthermore, the present transition pressure of 15 GPa agreed with experimental report. 7 ) 
Almost all MD simulation that graphite changes to diamond at comparable pressure to the 
experiment by the other researchers employed quantum chemical theory to calculate energy 
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and force of atoms. 14 ~ 22 ) It's interesting to reflect that the present new potential is classified 
into the category of classical mechanics though it can demonstrate the transition of graphite 
to diamond at proper pressure. 

5. Conclusions 

In the present work, we were aware that the region why sp 3 rich amorphous carbon ma- 
terials are not well generated by MD simulation is that the REBO2002 potential, which is 
major model, cannot give correct interaction energy in terms of amorphous carbon structures. 
To solve this problem, we developed new potential model, which was modeled according to 
that natural consideration that interaction among atoms consists of attraction and repulsion 
between nuclei and electrons. Then, electron distribution was simplified as a point charge 
located on places between atoms and same place to nuclei. Consequently, using the present 
new potential, we could generate sp 3 rich amorphous carbon deposits. In addition, we could 
demonstrate the transition from graphite to diamond at the pressure of 15 GPa, which was 
corresponding to experimental report. 7 ) 

To advance the research for amorphous carbon deposition, the treatment of hydrogen 
atoms is necessary. As additive substance, catalyst or impurity, more atomic element is needed. 
Therefore, we should tackle to extend potential model aiming to treat all atomic element. 
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Figure Captions 

Fig. 1. Change of pressure in the MD simulation of transition from graphite to diamond. The 
solid line indicate the setting pressure in the Andersen's method. The long dashed, short 
dashed and dashed-dotted lines indicate pressure acting to material for each x, y, z di- 
rections. Their plotted values are mean values for 1.0 x 10 -13 s. 

Fig. 2. (Color online) Animation snapshots in the MD simulation for transition from graphite 
to diamond at pressure of 15 GPa. The green and red spheres indicate sp 2 and sp 3 
carbon atoms. The white lines are the side of simulation box following periodic boundary 
condition. 

Fig. 3. The sp 2 and sp 3 ratios, which is defined as the ratios of the numbers of sp 2 and sp 3 
carbon atoms to that of all carbon atoms, as functions of the elapsed time. 

Fig. 4. The sp 2 and sp 3 ratio of injected carbon atoms in deposits. The white and black rect- 
angles indicate sp 2 and sp 3 ratios evaluated by the MD simulation with the present new 
potential, respectively. The white and black spheres indicate these with the REBO2002 
potential, respectively. 

Fig. 5. The penetration depth of carbon atoms as a function of injection energy. The penetra- 
tion depth is mean of that for all injected carbon atoms except for sputtered one. The 
circle and square points indicate the calculated depths by using REBO2002 potential and 
the present new potential, respectively. 

Fig. 6. (Color online) The thickness of deposit as a function of the elapsed time. 

Fig. 7. The inter-atomic interaction energies calculated by the present new potential, the 
REBO2002 potential and DFT using the OpenMX in terms of atomic geometries at 
the moment of transition process from graphite to diamond in the MD simulation with 
the present new potential. 
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